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We performed variational calculation based on the multi-scale entanglemnt renormalization 
ansatz, for the antiferromagnetic Heisenberg model on a Shastry Sutherland lattice (SSL). Our 
results show that at coupling ratio J' / J = 0.687(3), the system undergoes a quantum phase transi- 
tion from the orthogonal dimer order to the plaquette valence bond solid phase, which then transits 
into the antiferromagnetic order above J^ J = 0.75. In the presence of an external magnetic field, 
our calculations show clear evidences of various magnetic plateaux in systems with different coupling 
ratios range from 0.5 to 0.69. Our calculations are not limited to the small coupling ratio region, 
and we are able to show strong evidence of the presence of several supersolid phases, including ones 
above 1/2 and 1/3 plateaux. Such supersolid phases, which feature the coexistence of compressible 
superfluidity and crystalline long range order in triplet excitations, emerge at relatively large cou- 
pling ratio {J' / J > 0.5). A schematic phase diagram of the SSL model in the presence of magnetic 
field is provided. 

PACS numbers: 75.10.Jm, 75.10.Nr, 75.40.Mg, 75.40. Cx 



I. INTRODUCTION 

Geometrically frustrated interactions and quantum ef- 
fect often lead to exotic orders such as, quantum spin 
ice on a pyrochlore lattice spin liquid state on Kagome 
lattice^ and fractionalized magnetization plateaux^ be- 
cause strong fluctuations suppress stabilization of classi- 
cal or trivial magnetic orderings. A Shastry-Sutherland 
lattice (SSL)^ is one of such frustrated systems and the 
Hamiltonian of the SSL model can be written as 

H = J^S^-Sj^J' S,-S,-/i^^f, (1) 

where (i, j) and {{i^j)) stand for intra-dimer and inter- 
dimer site-pairs on the SSL as shown in Fig. (TJ panel 
(d). SrCu2 (303)2 has provided a good stage for study- 
ing the magnetic properties of the SSL model.^ The crys- 
tal structure of this compound is characterized by a lay- 
ered structure of CuBOs and Cu^+ ions carry S = 1/2 
spins 4 From magnetic susceptibility measurements^ it 
has been reported that the ground state is spin gapped 
and the coupling ratio is estimated as / J ~ 0.635.^-^ 
This compound shows anomalous properties in magnetic 
fields owing to the magnetic frustration; the magnetiza- 
tion curve exhibits a multi-steps structure. In addition to 
the plateaux with magnetization at 1/8, 1/4 and 1/3 of 
the saturation value, recent torque experiments sug- 
gested existence of 1/6 and other plateaux with more 
exotic fraction numbers ^i^i^^ 

Quite recently, rare-earth compound a^^i^^ — with the 
SSL structure were also studied in experiments, showing 
exciting magnetic properties may or may not differ from 
the original material SrCu2 (303)2 • While the inter- 



dimer and intra-dimer interaction strength differs from 
material to material and long range interactions can be 
crucial to understand magnetic properties in the rare- 
earth compound cases, the SSL model with frustrated 
Heisenberg interactions is still the key point to under- 
stand the physics of magnetic properties of these mate- 
rials. 

Theoretically, the ground-state phase diagram of the 
Hamiltonian ([1]) has also been in debate for this decade. 
The remarkable feature is that, in the limit jy J 0, the 
direct product state of singlets on J' bonds becomes the 
eigen state and the energy of this orthogonal dimer (OD) 
state can be given exactly by £^ = —3/8 J. Oppositely, 
the system becomes equivalent to the antiferromagnetic 
(AF) Heisenberg model on a square lattice in the limit 
J' / J 00. Thus the conventional AF state is realized 
with a gapless dispersion. This fact leads to naive expec- 
tation that a quantum phase transition between these 
two phases take place when J' / J is changed. Although 
the quantum phase transition including the possibility of 
intermediate phase has been intensively studied by an- 
alytical and numerical approaches ^li^—i^ conclusive 
result has not been established yet. 

In fact, obtained results strongly depend on ap- 
proaches. More specifically, the high-temperature- 
series expansion, early exact diagonalization and 
PEPSiiii^ii^i^ suggested the direct transition from the 
OD phase to the AF phase. In contrast, Schwinger 
boson meanfield theory suggested the existence of he- 
lical ordered state as an intermediate phase. Koga and 
Kawakami predicted that the other spin gapped phase 
exists for 0.677 < J' / J < 0.86.^^ and it can be ex- 
pressed by the alignment of resonating singlet dimers on 
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the plaquette without J bond, namely plaquette singlets. 
This prediction was supported by the other numerical 
exact diagonalization study on up-to 32-site systems 
More fascinating possibility was predicted from the re- 
sults of dimer-field theory in ref. |20|; weakly incommen- 
surate spin-density wave state is stabilized in the inter- 
mediate phase. The modified cluster expansion with self- 
consistent perturbation by Hajj and Malrieu suggested 
that the lowest energy state is either plaquette singlet or 
columnar dimer states The reason of the above long 
historical discussions on the intermediate region is caused 
by the fact that quantum fluctuation leads to numer- 
ous candidate states, all of which have similar energies. 
Therefore one may hope unbiased calculations such as 
quantum Monte Carlo (QMC). Since QMC calculation, 
which is very helpful to settle problems if possible, suffers 
from the negative sign problem in the present case, this 
approach will be excluded. 

Similar situation is also seen in magnetic properties 
when external field h exists: multi-plateau physics. Ex- 
act diagonalizations^^ and other pioneer studies^^ have 
established the existence of 1/2, 1/3 and 1/4 plateaux, 
while details of the small magnetization region were not 
very clear. In a recent paper by Sebastian et al.^ the 
possibility of fractional magnetic plateaux at 1/q {q = 
integer) was discussed in analogy to the quantum Hall 
effect. Several authors predicted many possible plateaux 
with further exotic fraction numbers like such as 1/9 
and 2/15^^ They focused on the strong dimer limit 
J' <^ J case, and inter-dimer interactions J' were treated 
as perturbations ^^JiS.^ In such case, triplet excitations in- 
duced by external magnetic fields can be regarded as 
hardcore bosons and the interactions on plaquette with- 
out J bonds can be considered as a combination of repul- 
sion interaction and hopping process. Although many 
believes that these plateaux shall persist in the large cou- 
pling ratio {J' / J > 0.5) region, up to best of our knowl- 
edge, no unbiased numerical method has thoroughly con- 
firmed the result in such situations. 

What is more interesting is the proposal of possible 
supersolid phases which may appear in the presence of 
large external field^^ The supersolid phases have been 
confirmed in various ID and 2D frustrated models, 
and it has been well established that correlated hop- 
ping (bosons hop when another boson is present as its 
neighbor) plays an important role,^^ since single boson 
hopping is strongly suppressed by frustration in the SSL 
model. Due to the strong competition between repul- 
sion and hopping terms, super soHd phases are unlikely 
to emerge as the ground state at small coupling ratio re- 
gion (jyj < 0.5). As a result, super solidty has not been 
well investigated in early calculations which treat J' as 
simple perturb at ion 

In this paper, we study the magnetic properties of the 
SSL model, using variational method based on the en- 
tanglement renormailzation (ER) ansatz. This method 
is more popularly known as "multi-scale entanglement 
renormailzation ansatz" , or simply MERA, which is pro- 



posed by Vidal a few years ago.^^ It is a promising numer- 
ical method based on merging the idea of tensor product 
state ansatz and renormailzation group theory. The key 
feature of MERA is the introduction of so called "disen- 
tanglers" which targets on reducing local correlations, so 
that long range entanglement can be captures with lim- 
ited computational cost at a coarse grained level. In our 
ER based calculation however, we simply use one level 
of coarse graining, instead of multiple levels in MERA 
calculations. We employ multiple neighboring unit cells 
in neighborhood of the unitcell we calculate, to study 
thermodynamic properties with meanfield approxima- 
tion. This method is also proposed by Evenbly and Vidal 
under the name "finite range MERA" (FRMERA),^^ due 
to the fact that only entanglement within finite range up 
to the size of each unticell is well captured. We will dis- 
cuss details of FRMERA and justify the reason we use 
this method in our ER calculation. 

We want to highlight that unlike previous studies men- 
tioned above, our calculation does not enforce any bias to 
the original Hamiltonian. The ground state is obtained 
from pure variational calculations, as the state with low- 
est energy for the input Hamiltonian (i.e., we provide 
upper boundary for the ground state energy.) Our cal- 
culation focuses on the relatively large coupling ratio re- 
gion, which has not been well investigated, since we are 
not limited to the small coupling ratio limit J' <C J. 

In the absence of magnetic field, we find strong evi- 
dence of the existence of the plaquette valence bond solid 
(VBS) phase in the intermediate coupling ratio region 
above the first phase transition point J' / J = 0.687(3). 
When the external field is turned on, we confirm the all 
the major (well accepted) plateaux such as 1/4, 1/3 and 
1 / 2 , as well as minor ones including 1 /6 , 1/8 and 1/9. We 
also observe a strong evidence of several supersolid orders 
above 1/3 and 1/2 plateaux, as proposed by Momoi et 
al.^^^ which has not been discussed in previous numeri- 
cal studies to best of our knowledge. Our results show 
that the 1/3 plateau only appears above J' / J = 0.5, and 
the 1/2 plateau and its corresponding supersolid phase 
persist until the largest coupling ratio J' / J = 0.687 be- 
fore the phase transition to the plaquette order. Our 
schematic phase diagram is shown in Fig. [6l Sec. IIVI 

The remaining parts of this paper is organized as fol- 
lowing: In Sec. ini we discuss the ER method and tensor 
network structure we designed for the SSL model. In 
Sec. results of the ground state in absence of mag- 
netic field are presented. We show magnetic behaviors of 
the SSL model from our calculation, including plateaux 
and supersoHd in Sec. |IVl Discussion and summary are 
presented at last. 



II. MERA METHOD AND TENSOR NETWORK 
STRUCTURE 

The structure of tensor network plays an important 
role in ER calculations. Although there are no strict 
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FIG. 1: Illustration of tensor netwrok structure used in ER 
calculation. 

Panel (a): The shastry sutherand model plotted on a square 
lattice. Intra-dimer (J) and inter-dimer (J') couplings are 
shown as solid and dashed line, respectively. In the entan- 
glement renormalization procedure, 8 original spins included 
in each domain (shown as shadow region) are coarse grained 
into one (shown as blue dot) by a blue isometry tensor. Two 
sets of disentanglers (red and green squares) are applied on 
boundaries between different domains. 

Panel (b) : two sets of disentangler tensors and one set of isom- 
etry tensors are applied in a consecutive manner, following the 
coarse graining direction from bottom (starting hamiltonian) 
to top (coarse grained system). Lower and upper legs stand 
for in and out indices, respectively. 

Panel (c): a final density matrix is applied on the coarse 
grained system (indicated by blue circles). For an unitcell 
(denoted by shadows in the background) with 4x6 spins, 3 
coarse grained sites are involved in the density matrix, which 
is indicated by a light triangle. For unitcell shaped 6x8, 
6 coasre grained sites are inculded, by doubling the length 
along the y direction. In FRMERA calculation, neighboring 
unitcells (shadow regions) are used as meanfield environment 
to simulate system in the thermodynamic limit L ^ oo. 
Panel (d): Unitcells consist of 4 x 8 spins, or equivalently 
4 coarse grained sites. This unitcell is specially chosed to 
perserve the X-Y symmetry of the square lattice. 



rules for constructing ER tensor network for a specific 
model, one guideline shall be taken into consideration. 
The advantage of ER based method comes from the dis- 
entangler tensors, which reduce local correlations so that 
the density matrix of the coarse grained (CG) system 
is capable of capturing long range entanglement of the 
original spin model with relatively small virtual dimen- 
sions X- Naturally, all interactions across boundary of 
the coarse graining domain should be treated by disen- 
tangler tensors. Relatively speaking, in most frustrated 
two dimensional quantum spin systems, it is common to 
have more interactions involved with each site, compared 
to non- frustration spin models. As a result it is impor- 
tant to design an ER tensor network suitable for the SSL 
model that we are interested in. 

In Fig. [H we show an illustration of the ER tensor 
network structure we chose for this calculation. Notice 
that the network is constructed on the square lattice's 
geometry, instead of the SSL's. We believe that it does 
not affect the result much in the small coupling ratio 
limit, as the ground state can be considers as decoupled 
dimers separated on the lattice. Moreover, such a choice 
may aid us in the large coupling ratio region where we 
are most interested in. This geometry allows us to easily 
utilize our tensor structure to study other frustrated lat- 
tice which can be mapped onto square lattice (i.e., the 
checker-board model). 

Overall, 8 sites on the original square lattice are coarse 
grained into one site on the renormalized lattice, which 
is another square lattice tilted by 7r/4. Coarse-graining 
sites obtained from isometry tensors and the correspond- 
ing domain containing the 8 original spins are shown as 
blue circles and shadowed regions, respectively, in Fig. [TJ 
panel (a). Two layers of 4-sites disentanglers (colored by 
red and green) are applied consecutively on the boundary 
of coarse-graining domains, before the isometry tensors 
are applied. Overlook the tensor network on top of the 
SSL, it is clear that all intra and inter dimer interactions 
are covered by at least one tensor of the network. As a re- 
sult, all interactions across the boundary of CG domains 
are treated by disentanglers. 

The coarse grained lattice is then described by an top 
density matrix. Obviously, the dimension of the coarse- 
grained spin X cind the number of CG spins covered by 
the density matrix determine computational cost of our 
calculation. In practice , we are able to treat up to 6 
coarse grained sites (4x8 spins) per unitcell. 

As mentioned previously, unlike the MERA method, 
we only apply one level of renormalization to the SSL in 
our calculation. Instead of limiting ourselves to finite-size 
systems with periodic boundary conditions, we choose 
FRMERA to construct multiple identical neighboring 
unitcells. In Panel (c) of Fig. [H we illustrate an example 
when the unitcell contains 4 x 6 = 24 spins (shadowed 
region in the background). In such a case, the density 
matrix of each unitcell (represented by light triangles in 
Fig. d]) is not correlated to its neighbors at all. As a 
result, the FRMERA method is obviously related to the 
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meanfield approximation in the way how we simulate sys- 
tem in thermodynamic hmit. 

The reason we choose FRMERA instead of the ordi- 
nary MERA is twofold. First, although applying mul- 
tiple levels of renormailzation can significantly increase 
the system size we can study in MERA, associated cost 
is very high. This is due to the fact that sufficiently 
large virtual dimensions x is required for each site in 
the coarse grained system in order to capture sufficient 
long range entanglement in higher level of renormailza- 
tion. Moreover, we are not specially focusing on quantum 
critical region whereas scale-invariant physics play an im- 
portant role (where Multi-level renormailzation is impor- 
tant and meaningful). Rather, we pay special attentions 
to the ground state of SSL model with and without pres- 
ence of external magnetic field. In such a case, as long 
as the unitcell is large enough to accommodate candi- 
date ground states and its associated entanglement, we 
can simply use ER with meanfield approximation which 
is significantly cheaper in computational cost compared 
with MERA. 

As a matter of fact, FRMERA turns out to be bene- 
ficial in studying the magnetic properties of SSL model. 
The SSL model shows intriguing magnetic properties in 
response to external magnetic field h. Besides 1/3,1/4 
and 1/8 magnetic plateaux observed in early experi- 
ments, there are theoretical studies as well as possi- 
ble signature from experiment supporting existence from 
plateaux with small fraction numbers, such as 1/6, 1/9 
and 2/15^^2*21 One difficulty to study these magnetic 
plateaux computationally lies in the fact that they only 
appear as ground states if the unitcell itself is compat- 
ible with the fraction number. For example, plateaux 
with fraction number 1/3, 1/6 and so on can not be well 
described in ER calculations performed on 4 x 8 unitcell, 
nor can it be even with multiple renormailzation applied 
(due to periodicity). Obviously, it is increasingly difficult 
to capture plateaux with even smaller fraction numbers, 
which require larger unitcells to settle in. The FRMERA 
method offers a semi-solution to this problem. Namely, 
one can further sort unitcells into different categories, 
each contains a different set of tensors. Following this, 
we can effectively create a larger super-unitcell, which 
combines one unitcell per category, at a computational 
cost linear to the size (number of unitcells included), we 
show such super-unitcells as background in Fig. 21 Notice 
that unitcells shown in different colors belong to different 
categories. As one example, in the last panel, the super- 
unitcell is considered as a combination of 3 unitcells, one 
per each color, so that the ER calculation is capable to 
capture the 1/9 plateau. 

It is worth mentioning that unitcells included in the 
super-unitcell do not have correlations with each other at 
top level, since it is just a product state like a meanfield 
approximation. On the other hand, correlations can still 
pass through in lower level tensor network, so that we are 
still able to study plateaux with small fraction numbers 
which would be impossible to find in ordinary unitcell. 
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FIG. 2: Upper panel: AF and dimer order parameters, M and 
D, as a function of coupling ratio J. Results obtained from 
unitcell sized 4x6 and 4x8 are shown in black and red (light) 
colors, respectively. A strong first order transition occur at 
J' / J = 0.687(3), where the system enter the plaquette VBS 
order which features finite dimer order parameter D. The 
transition point between the plaquette order and the standard 
AF order is not clear, as the crossing point between M and D 
shifts significantly as the unitcell size is increased. Energies of 
ground state (normalized by J') are shown in the inset, with 
dashed lines plotted as "guide for eyes" for phase boundaries. 
Energies obtained from larger unitcell 4 x 8 is slightly better 
compared with those obtained for 4x6 unitcell. Lower panel: 
local bond energy profile for the OD order and the VBS order. 
Thickness of the bond corresponds to the correlation between 
two spins {Si ■ Sj). 



or naively increase number of renormailzation steps. 

The dimension of coarse grained virtual spin, ef- 
fects our calculation results. For unitcells consists of 3, 
4 and 6 coarse grained spins per density matrix, (which 
correspond to 24, 32, 48 original spins,) we are able to 
perform calculations up to x = 9? 9 and 4, respectively. 
It is necessary and useful to check the convergence (or 
evolution trend) of physical quantities as a function of 
the dimension x- We will discuss behaviors of evolution 
trend when needed in following discussions. 



III. GROUND STATE OF SSL MODEL IN 
ABSENCE OF EXTERNAL FIELD 

Without external magnetic field, the true ground state 
of the Heisenberg model on the SSL remains to be contro- 
versial, but there are two trivial reference limits for the 
ground state. In the limit J' <C J, the system can simply 



5 



be considered as a product of decoupled singlet dimers, or 
commonly referred to as the OD model. In another limit, 
the system reverts back to the non-frustrated Heisenberg 
model on the square lattice, with the AF order being the 
ground state. The debate is mainly focused on the inter- 
mediate region where quantum fluctuations are large and 
sign-problem due to frustration of SSL model prohibits 
effective numerical methods for quantum magnets to be 
performed. In Fig. O we show energy (inset) and order 
parameters as a function of coupling ratio J' / J . 

To discuss the possibility of intermediate phase, 
namely plaquette phase, we introduce two estimators as 
the order parameters. The ordinary AF order parameter 
M is defined as 

M = ^Ml + Ml + Ml, (2) 

i 

(4) 

and M^, are defined likewise. The order parameter 
D for the valence bond solid order is defined as 

= ^E(s.-Si+.)(-ir% (5) 

i 

= ^E(Si-s.+s)(-ir% (6) 

i 

D = {Dl + Dl)i, (7) 

where x^, yi are coordinates of site i, and (z, j) stands for 
nearest neighbor pairs on the square lattice. and Dy 
are not direct representations of plaquette state, but they 
are ideal estimators because we can distinguish colum- 
nar and plaquette VBS order by evaluating the differ- 
ence between x and y components of the dimer order, 
D/\ = \ \Dx \ — \Dy\\^ which is also shown in Fig. [2j 

It is clear that the system is in the OD phase in the 
small couphng ratio region, because the energy can be 
given exactly by E = —3J/S. In fact, AF and dimer 
order parameters are strictly zero, and the energy of the 
system is linear to J^ 

As J' / J increases, the ground state energy show a cusp 
at jy J = 0.687, and the dimer order parameter D jumps 
to a finite value keeping the value D/\ ~ 0. (see red 
squares and diamonds in Fig.[2j) The system clearly un- 
dergoes a first order phase transition from the OD phase 
into the plaquette VBS phase. The estimated transition 
point is in good agreement with that in previous works. 
As further evidence of the plaquette VBS order, we eval- 
uate the local bonding energy and show the pattern in 
the lower panels of Fig. [2l The bond energy profile also 
strongly indicates that the ground state in the interme- 
diate region does not break 7r/2 rotational symmetry. 

Results of FRMERA calculation are dependent to the 
size of unitcell we use. In Fig. [21 we show two sets of data 
which correspond to calculations based on two types of 



unitcells sized 4x6 and 4x8. The dimer order pa- 
rameter D increases when the unitcell is enlarged. The 
4x6 unitcell clearly breaks the rotational symmetry of 
the original model, whereas the fan-shaped 4x8 unitcell 
does not (both shown in Fig. [T]). As a result, in the in- 
termediate region just above J' / J > 0.687, D/\ remains 
small but finite in 4 x 6 unitcells' calculations, but turns 
out to be in the 4x8 sized unitcells' case. We conclude 
that plaquette VBS state, which preserve 7r/2 rotational 
symmetry of the square lattice, emerges as the ground 
state when the system undergoes a phase transition out 
of the OD order. 

In contrast to the phase transition point between the 
OD phase and plaquette VBS phase, the transition be- 
tween the AF order and plaquette VBS phase seems not 
to be clear due to finite size effect of our unitcell. The 
value of D drastically decreases around J ^ 0.75 and 
M starts to develops inversely. However, the crossing 
point of two order parameters is shifting to higher J' / J 
value, when the size of the unitcell is increased. In fact, 
order parameters are slightly affected by dimension x 
(i.e., not fully converged) for 0.85 > J' / J > 0.75. More- 
over, both order parameters remain finite in a relatively 
large region near the transition point. As a result, solely 
based on available result, we can not clearly determine 
the nature of the phase transition between the plaque- 
tte VBS order and the AF order. Further calculation 
is required for precise discussion of the phase transition 
point. 

It is necessary to discuss the validity of the ground 
state we obtained from ER calculation. First, it is natu- 
ral to expect that MERA and ER methods are biased in 
the sense that the tensor network (including the way we 
choose coarse grained lattice) has a spacial structure. Ac- 
tually, the ER (and MERA) procedure probably breaks 
spacial (translational and rotational) symmetries inher- 
ing in the original Hamiltonian. As a result, optimized 
wave functions of the ground state may contain unwanted 
effective perturbation introduced from the tensor struc- 
ture. Naturally, if the dimension x of coarse grained 
system is sufficiently large, the effect of such bias can 
be neglected, and the solution of the ER calculation be- 
comes exact. In practice, x used in our calculation is 
limited, to avoid explosion of simulation cost. Such lim- 
ited X will cause severe problem to determine the nature 
of the ground state in the critical region, or to locate the 
phase transition point, where the gap between energies 
of different candidate orders are closing. 

What we can do is to reduce such bias by choosing an 
ER tensor structure that is compatible to most candi- 
dates of ground states. For the present case, we specif- 
ically chosen a tensor network structure with even peri- 
odicity (i.e., the coarse graining bulk contains 8 sites) as 
discussed in the previous section, because the plaquette 
VBS state (4 site per unit cell) is one of the candidate 
ground states at a zero field. Apparently, the network 
adopted has suitable structure for the wave function of 
most candidates, namely the AF state, the OD state, and 
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the VBS state. Unfortunately the ER structure we use 
has difficulties dealing with incommensurate spin den- 
sity wave (SDW) state, as common for most numerical 
method based on a finite unitcell with periodic bound- 
ary conditions. If we can apply mutiple- level ERs on the 
SSL, the incommensurate SDW state may be captured 
approximately as in ref. |40|. However, due to very high 
computational cost, it is difficult in our ER. To summa- 
rize, we believe that our result is not much affected by 
the bias introduced from tensor network structure away 
from the critical region, although we can not rule out the 
possibility of the incommensurate SDW phase. 

Secondly, the ER calculation we performed is essen- 
tially a variational approach. A well known feature of 
variational calculation is that optimization may converge 
into local minimums, instead of the global minimum - the 
true ground state. In terms of energy landscape, the OD 
phase and the plaquette phase (which is fourfold degen- 
erate on the square lattice) correspond to far separated, 
sharp and stable local minimums. In other word, near the 
critical region, optimization calculations may end in both 
candidate states with finite possibilities, and it is impos- 
sible to "tunnel" between two different phases whence 
the optimization has reached the local minimum. Hence, 
this is also a strong indication that the phase transition 
between the OD order and the plaquette VBS order is 
of ffist order. A jump in order parameters and a sudden 
change of slope in energy can be observed in Fig. [2j In 
contrast, in the critical region between the VBS and AF 
phase, the two local minimums corresponding to the re- 
spective phases are shallow, not well separated and show 
tendency of merging with each other at the critical point. 
Therefore the slope of energy does not change drastically, 
while D and M evolve smoothly across the phase bound- 
ary. This is an indication of a second order phase transi- 
tion, but we cannot rule out the possibility of cross over 
because of the severe finite-size dependence. 



IV. MAGNETIC PROPERTIES OF SSL MODEL 

What makes the SSL model and its realization 
SrCu2(B03)2 fascinating is its magnetic property in a 
magnetic field. It is well known from experiments that 
the coupling ratio is J' / J ~ 0.65, where the ground state 
is given to be OD state (product of singlets). In the 
presence of an external magnetic field, polarized triplets 
with Sz = I are created, which can be treated as hard- 
core bosons moving on the effective square lattice. Due 
to strong frustration in the lattice, single triple hopping 
is strongly suppressed, and effective repulsion between 
triplets governs the physics in the system when the cou- 
pling ratio satisfies jyj < 0.5. Naturally, magnetic ex- 
citations form long range order and then, the system is 
in a magnetic plateau: an incompressible Mott phase. 

Pioneer theoretical works conffimed several magnetic 
plateaux. Among them, 1/3 and 1/2 are most well ac- 
cepted, while plateaux in lower magnetization are still 




FIG. 3: Magnetization of the system in response to exter- 
nal field h. Results obtained at three different coupling ra- 
tio f/J = 0.5,0.635,0.667 are shown. 1/2, 1/3 and 1/4 
plateaux are clearly shown. 1/9, 1/8 and 1/6 plateaux are 
also observed. Open symbols denote volume averaged off- 
diagonal spin component at three coupling ratios. 
They are indications to the supersolid phase. Inset: The ener- 
gies of variational local minimum branches that corresponds 
to states with 1/2, 1/3 and 1/4 plateaus' dimer structures. 
Cross points between two curves stand for jumps of magneti- 
zation. 



less clear. Recently, Sebastian et.al}^ have suggested 
1/9, 1/6 plateaux and the other exotic ones at frac- 
tionalized magnetization from the quantum Hall anal- 
ogy. Other theoretical woik.^^^-^-^ also successfully 
confirmed the existence of these plateaux based on per- 
turbative approaches from small J' / J limit. The result of 
our calculation is largely consistent with previous studies, 
but we have to highlight the following point: our calcu- 
lation is not restricted to the small J' / J limit, since we 
do not treat J' interaction as a perturbation of the sin- 
glet product state. As a matter of fact, we simply input 
our hamiltonian in the original square lattice, as shown 
in Fig. [H and the magnetic plateaux are obtained as the 
ground state of our pure variational calculation. As a 
result, besides small J' / J ratio (< 0.5), we also investi- 
gated large coupling ratios close to the phase transition 
point (into the plaquette VBS order). 



A. Magnetic plateaux 

It is clear judging from Fig. [3l that we obtain well ac- 
cepted 1/4, 1/3 and 1/2 plateaux in the magnetic curve 
when external field is applied. In addition to these three 
plateaus, minor plateaux such as 1/6, 1/8 and 1/9 are ob- 
served as well In Fig. 21 we plotted correlations between 
neighboring spin pairs (•) (for both J and J' bonds), 
as well as {Sz) for all local spins, for minor plateaux 
(similar information for 1/2 and 1/3 can be read in 
Fig. [5]). These dimer patterns are consistent with previ- 
ous studies Notice that color and width of each 
dimer we draw stand for spin correlations between two 
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sites forming the dimer, where red and green stand for 
ferromagnetic (triplet) and AF (singlet) correlation, re- 
spectively. The magnitude of Sz for each local spin is 
also plotted as the radius of circle on each site. 

It is worth mentioning that the 1/8 plateau appears 
only around J' / J = 0.5 and the local moment distribu- 
tion (S^) is of the so called rhomboid type, although the 
energy difference between rhomboid type and square type 
is extremely small. For jyj > 0.5, we find that the 1/8 
plateau state (rhomboid type configuration) has higher 
energy compared with other plateaux. This discrepancy 
may arise from further long-range interactions. 

Another interesting point is that the 2/9 plateau state, 
proposed in ref. HH, which can be covered by our super- 
unitcell tensor network wave function (and emerge as lo- 
cal minima in our variational calculation), turns out to 
be worse in energy compared to other low magnetiza- 
tion plateau states (1/4,1/6) in relevant parameter re- 
gion. This fact indicates that the PCUT treatmenli2S 
seems to be already broken down for jyj > 0.5. 

Again here we discuss the validity of our variational 
ground states. It is obvious that plateau states with dif- 
ferent fraction number correspond to local minimums, in 
terms of energy landscape in our calculation. In the pres- 
ence of relatively large external field, (i.e., m/rris > 1/3, 
while rris in the saturation magnetization), we find there 
are two major groups of local minimums: one corre- 
sponds 1/2 plateau and the other is 1/3. These two lo- 
cal minimums are stable, sharp and far separated, which 
means tunneling from one plateau state to another is 
nearly impossible to happen, especially at the later part 
of the optimization. The reason is simply that these two 
plateau states are characterized by different periodicity. 
We use the existence of stable local minimums to our ad- 
vantage. We can obtain energy of a certain plateau state 
as a monotonic function of field /i, if we choose to load 
tensors from previous calculations (that is know to host 
the plateau state) instead of starting from random initial 
condition. Phase boundaries are then determined by the 
crossing point of energy curves of different plateau states, 
as shown in inset Fig. [3l Each crossing point corresponds 
to a jump in the magnetic curve, which is a strong first- 
order phase transition. Similar techniques were also used 
in other variational method as well^^ However, in lower 
magnetizations (m < 1/4), the energy gap between dif- 
ferent plateau states are fairly close and the phase bound- 
aries between different plateau states are less clear, be- 
cause the energy for each plateau state is sensitive to 
changing of unit cell geometry and other factors. It is 
also necessary to point out that the dimension x of 
coarse grained system in our calculation has very limited 
impact on results of magnetic plateaux. In fact, all mea- 
sured observables are not much dependent on x, and the 
magnetization tends to reach the correct fraction number 
even at smallest dimension x = 2. 
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FIG. 4: Two sites spin-spin correlations i^Si • Sj) for J and J' 
bonds plotted on the square lattice. From red color to green, 
the correlation turns from ferromagnetic (i.e., triplets) into 
AF (singlets) . The width of each plotted bond is proportional 
to the magnitude of correlation, (i.e, very thin line and dim 
color represents near-0 correlations.) Open/closed circle at 
each site represents polarized/anti-parallel local spin, where 
radius of each circle stands for magnitude of {Sz)- Shadow 
region in the background shows the arrangement of unit cells. 
Top two panels: correlation patterns obtained in 1/4 and 1/8 
plateau states. Bottom two panels: patterns obtained in 1/6 
and 1/9 plateau states. Large super-unitcells, which consist 
of one unitcell per category (represented by different color), 
are used to accommodate plateau states with small fraction 
numbers. 



B. Supersolid states 

In addition to multi plateaux, it was predicted that 
spin supersolid phases appear in the SSL model when 
magnetic fields exists ^^i^^ It has been pointed out that 
correlated hopping of triplets is crucial to stabilize the 
supersolid phases. The correlated hopping (i.e., triplets, 
which are regarded as bosons, hop from one dimer to an- 
other when a neighboring boson is present) contributes 
as the main source of kinetic energy in the system, in con- 
trast to single triplet hopping which is largely suppressed 
by frustration. Such correlated hopping can lead to su- 
persolid phase, where diagonal (crystalline solid) and off- 
diagonal (superfluidity) long range order coexist For 
the SSL model with considerably large coupling ratio 
jyj, supersolid phases are predicted to appear above 
1/3 and 1/2 plateaus in the high fleld region. 

Quite recently, high fleld measurements on 
SrCu2 (303)2 over 100 [T] have been achieved and 
the presence of the 1/2 plateau has been conflrmedj^ 
In real compound, the Dzyaloshinski-Moriya (DM) 
interaction, which breaks U(l) symmetry extrinsically, 
can not be ignored, therefore superfluid transition is 
not prohibited in the compound, when the temperature 
decreases. Fortunately, the magnitude of the DM inter- 
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FIG. 5: Dimer patterns obtained for 1/2, 1/3 (top two pan- 
els) and P3 (bottom two panels) supersolid states. The width 
of each dimer is proportional to the magnitude of correlations 
between two spins. From AF to ferromagnetic correlations, 
the color of each dimer can change from green to red con- 
tinuously. Arrow on the each dimer stands for the in-plane 
component of the coherent vector ($,G), and the radius of 
circle on each dimer is proportional to the vector component 
normal to the plane. 



action in SrCu2(B03)2 was estimated as I^dm ^ 2[K}^ 
and it is about 0.2% of the strongest coupling J. This 
allows us to expect that the weak DM interaction would 
fix the phase of off-diagonal components. Note that the 
transition related with superfuildity is not recovered 
due to the intrinsically symmetry of spin Hamiltonian. 
Beside this point, theoretically, the accuracy of previous 
studies' evaluation of correlated hopping effect in large 
jyj region lacks, and to the best of our knowledge, 
supersoHd phase are not clearly observed in unbiased 
numerical studies for the SSL model. 

From Fig. O it is clear that the magnetization m re- 
spond to the external field continuously (gapless mode) 
just above the 1/3 and 1/2 plateaux, strongly indi- 
cating existence of supersolid phases. To further sup- 
port this idea, we notice that all dimers in our sys- 
tem can be considered as a mixing of singlet state and 
triplet excitation. Using a classical coherent vector = 
cos^sinO, sin^sinQ^cos^), we can represent a dimer 
state as: 



|i.,e)=e'* sin( I 



|s) + e-f cos(|)|t), (8) 



where \t) and 1^) stand for triplet and singlet states, re- 
spectively. Only when superfluidity of triplet excitations 
exist in the system, O can take angles not equal to 0,7r. 
It is easy to show that under such circumstance, the off- 
diagonal components of original spins in the lattice take 



non zero values as: 



{St) 



2V2 



sine, {S+) 



2V2 



sinB. 



(9) 



As a result, off-diagonal elements are anti-parallel for the 
two sites 1 and 2 on a dimer bond. 

We measure off-diagonal elements of original spin 
as the order parameter of supersolid phase. Results are 
shown in Fig. [3l indicating presence of strong supersolid 
phase above 1/3 (for large coupling ratio J'/ J = 0.635 
and 0.667) and 1/2 (for all J'/ J parameters we tested) 
plateaux. Notice that the 1/3 supersolid phase appears 
only in large coupling ratio region, due to the fact that 
for large J^J^ the repulsion between dimers are largely 
reduced, while related correlated hopping remains finite. 
On the other hand, for relative small J' / J = 0.5, the 1/3 
supersolid does not appear as the ground state in our 
calculation. 

The geometry of crystalHne pattern of the triplets plays 
an important role in the correlated hopping process, as 
pointed out by several studies. ^^^^^^^^ More specifically, 
for both 1/3 and 1/2 plateau states, triplet excitations 
(and the remaining singlet dimers) form stripes in (tt, tt) 
direction of the SSL, which corresponds to horizontal 
and vertical lines in our illustrations (Fig. [5j). In the su- 
persoHd phase, additional triplets can hop both parallel 
to the stripes (i.e., flow inside the "canal" between two 
triplet stripes) and perpendicular to it ("jump" across 
the stripe). It has been pointed out that hopping matri- 
ces with regard to hopping inside and across a stripe of 
triplet are of opposite sign. As a result, the superfluid 
components (presented by angle ^ in ([8])) is parallel in- 
side stripes, and anti-parallel across stripes. 

Using the local density matrix obtained from MERA 
calculation, we are able to calculate the superfluid com- 
ponent for each local dimer, as 



Z = \t){t\-\s){s\, 
X = \t){s\ + \s){t\, 
Y = ii\t){s\ - \s){t\). 



(10) 

(11) 
(12) 



Resulted coherent vectors are plotted in Fig. O where 
arrows show the direction of superfluid in the XY plane. 
Our result is consistent with the expectation for both 
1/2 and 1/3 supersolid phases. Notice that although su- 
persoHd vectors form LRO, the direction of such ordered 
"magnetization" itself is not fixed, which can turn out to 
be any angle within the plane, since random tensors are 
used as the starting point of our variational calculation. 
Due to the ordering of superfluid component, the spatial 
symmetry in the supersolid phase is different compared 
with respective plateau order. One need to be cautious 
setting up the unitcell for numerical simulations of su- 
persoHd phases. 

Based on these results, we provide a schematic phase 
diagram in Fig.[6l with coupling ratio J' / J and external 
field h as parameters. At J' / J = 0.5, the phase sepa- 
ration point from our calculation is in good agreement 
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FIG. 6: We propose a tentative phase diagram focused on 
large coupling ratio region J' / J > 0.5. The 1/2 plateaux 
persist to the largest coupling ratio J' / J — 0.667 (close to 
the transition point to the plaquette VBS phase) we have ex- 
amined. 1/8 plateau only appear at J' / J = 0.5. Supersolid 
phase exist just above the 1/3 and 1/2 plateaux at large cou- 
pling ratio. Above 1/2 supersolid, another supersolid phase 
with period 3 emerge as ground state. 



magnetization keeps responding continuously to external 
field, and {S^) remains finite, as shown in Fig. [3l In 
contrast to early studies, this supersolid phase emerges 
in our calculation and features a crystalhne pattern of 
singlet-triplet description. As a matter of fact, it can 
be considered as a natural extension to the 1/3 plateau 
state, as the period remains to be 3 rows of dimers. More 
specifically, when field h is large enough, while one stripe 
of triplet remain unchanged, two remaining stripes of sin- 
glets start to covert to triplets, one row follows another, 
until all dimers are polarized. The conversion is in a 
continuous fashion, i.e., the correlations (SiSj) on the 
dimer change from ferromagnetic to AF smoothly. (See 
example in Fig. O) We refer to such supersolid phase as 
P3 supersolid. The reason that the P3 supersolid phase 
is favored in energy compared to 1/2 supersolid at large 
magnetization region is unclear. One possible reason is 
that during the smooth polarization from low magnetiza- 
tion to full saturation, (i.e., 1/3 — 2/3 — 1), the diagonal 
LRO is largely unaffected, which features a period of 3 
rows of dimers. Correlated hopping may also be easier 
in P3 crystalline patterns, which render P3 supersolid an 
advantage in kinetic energy. 



with the phase diagram from Capponi et.al.^ Notice 
that 2/9 plateau is not well established in our calcula- 
tion, as the energy of such LRO pattern is not favored 
in energy (although they do appear as local minimum in 
optimization). 

Since our method is not limited to small coupling ra- 
tio limit, we mainly probed large J' / J region. Our re- 
sult largely agree with the schematic phase diagram pro- 
posed in strong-coupling expansion calculation by Momoi 
et.alS^ with a few discrepancies. 

First of all, the external field region h where the 1/2 
plateau exists does decrease significantly when the cou- 
pling ratio J' / J is increased, but persist until the largest 
ratio we tested, J' / J = 0.687, which is the phase sepa- 
ration point to the OD order in the original SSL model 
with no external field. The main reason for the shrink- 
ing of 1/2 plateau has been revealed in previous studies: 
the strength of the nearest neighbor repulsion between 
dimers significantly reduces in large J^J region. It is 
however not clear whether the 1/2 plateau shall remain 
in the compound SrCu2(B03)2 (with jyj ^ 0.635), or 
materials with even larger coupling ratios. According to 
our ER calculation, the 1/2 plateau still gains its foot in 
the competition against supersolid phases at large cou- 
pling region. 

Secondly, the 1/3 supersolid phase does not appear 
at jyj = 0.5, rather, emerge at fairly large coupling 
ratio. This is slightly inconsistent with Momoi's predic- 
tion, but agrees with more carefully executed perturba- 
tion calculations 

Above the 1/2 supersolid phase, through a frist or- 
der phase transition, the system enters another super- 
solid phase (with different triplet pattern), in which its 



V. SUMMARY 

Using MERA method, we are able to calculate the 
ground state of the SSL model both in the absence and 
presence of external magnetic field. Our calculation 
does not make any approximation of the Hamiltonian 
itself, unlike some early studies, and the ground state 
we obtained is the result of pure variational calculation. 
Thermodynamic limit has been reached by meanfield- 
like approximations (FRMERA, although there is sub- 
tle difference). We have showed clear evidence of a 
plaquette VBS order intermediating the AF and the 
decoupled dimer phase. We have found various mag- 
netic plateaux in the presence of external field, includ- 
ing 1/2,1/3,1/4,1/6,1/8,1/9. We have also confirmed 
the existence of supersolid phases above 1/3 and 1/2 
plateaux, where superfiuidity (off-diagonal LRO) and 
crystalline pattern of triplet excitations coexist. Above 
the 1/2 supersolid, another supersolid phase which fea- 
tures period 3 dimer lattice structure persists all the way 
to saturation. 

We have to point out that the FRMERA method faces 
difficulties in dealing with incommensurate phases, due 
to finite unitcell, and lacks in description of long range 
entanglement that is far beyond the size of the unitcell. 
Nevertheless, it is a very effective numerical method for 
studying gapped quantum spin systems with frustration. 
We have to point out here that DM interaction, which 
is believed to have profound impact in the material's 
magnetic properties, is regrettably not included in our 
discussion. Also effect of higher-order-neighbor interac- 
tions have been not included either^^. We notice that 
it is practical to include them in the calculation, if our 
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MERA tensor network is modified accordingly in order 
to restrain the cost. We choose to leave these possibilities 
to future studies. 
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